rm(list=ls())
setwd("/Users/yuehou/Dropbox/KingReplication/Ongoing Work")
 
# DV: Govem_yh2    IV: leftgs (lag_leftgs,lag_leftgs2)


#############################
#############################
#### Loop with Beck-Katz ####
#############################
#############################

# NOTE: You must input the models (below) manually 

# This is a hack of zelig functionality, that allows for correct
# standard errors and coefficients, but also allows the use of
# non-zelig models (plm), etc
# It also allows the user to modify the imputed dataset by adding lags, etc


load("imputations23.Rdata")
m <- 8

# Loads initial imputations into dataframe
for (j in 1:m){
	assign(paste("data_",j,sep=""), amelia.out9$imputation[[j]])
}


for (j in 1:m){
	data_temp <- get((paste("data_",j,sep="")))
	

	# Reframe as a panel dataset
	library(plm)
	p_newset <- pdata.frame(data_temp)
	
	p_newset2<-subset(data_temp,dhasminwage==1)  #just for minwage
	p_newset2<-pdata.frame(p_newset2)

	p_newset$log_oecd_gdp <- log(p_newset$oecd_gdp)
	p_newset$lag_oecd_gdp <- lag(p_newset$oecd_gdp)

# Create lag terms here
	p_newset$lag_ia5010 <- lag(p_newset$final_ia5010)
	p_newset$lag_govem_yh2<- lag(p_newset$govem_yh2)

	p_newset$lag_leftc<- lag(p_newset$leftc)
	p_newset$lag_leftc_2<- lag(p_newset$leftc,2)
	p_newset$lag_leftc_3<- lag(p_newset$leftc,3)
	p_newset$lag_leftc_4<- lag(p_newset$leftc,4)
	p_newset$lag_leftc_5<- lag(p_newset$leftc,5)

	p_newset$lag_leftgs   <- lag(p_newset$leftgs)
	p_newset$lag_leftgs_2 <- lag(p_newset$leftgs,2)
	p_newset$lag_leftgs_3 <- lag(p_newset$leftgs,3)
	p_newset$lag_leftgs_4 <- lag(p_newset$leftgs,4)
	p_newset$lag_leftgs_5 <- lag(p_newset$leftgs,5)

	p_newset$lag_arm_generos = lag(p_newset$arm_generos)
	p_newset$lag_minwag_jf = lag(p_newset$minwag_jf)
	
	

# Transform generosity measure

	temp1 <- p_newset$soc_exp1 <- (p_newset$arm_generos * p_newset$oecd_gdp) - (p_newset$lag_arm_generos*p_newset$lag_oecd_gdp)
	
	temp2 <- p_newset$lag_arm_generos * p_newset$lag_oecd_gdp
	
	temp3 <- p_newset$soc_exp <- 100 * (temp1 / temp2)

	# Oil covariate
	p_newset$oilshock <- 0
	p_newset$oilshock <- (factor(p_newset$year)==1973 | factor(p_newset$year) == 1972 | factor(p_newset$year) == 1971 | 	factor(p_newset$year) == 1970)



	
 ###############################################
   

results<- plm(model="within",data=p_newset,govem_yh2 ~ lag_govem_yh2 + vk_corp:leftgs + leftgs + vk_corp:lag_leftgs + lag_leftgs + lag_leftgs_2 + vk_corp:lag_leftgs_2 +vk_corp + oecd_gdpgr + arm_unemp + arm_open + arm_finop + oecd_debt + oilshock) # Two Lag


results2 <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))
 
 
   # marginal effects (govem and corporatism)
   hkvalue<-seq(1,5,1)
   
   beta_interaction<- results$coeff["leftgs"]+hkvalue*results$coeff["vk_corp:leftgs"]
   
   covariance1 <-vcovBK(results,cluster=c("group","time"))["leftgs","vk_corp:leftgs"]
   
   
   se_interaction <- sqrt(results2["leftgs",2]^2+(hkvalue^2)*results2["vk_corp:leftgs",2]^2+2*hkvalue*covariance1)



	results <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))

	# Output vectors
	assign(paste("resultsc_",j,sep=""),results[,1])
	assign(paste("resultss_",j,sep=""),results[,2])
	
	assign(paste("resultsinteraction_",j,sep=""),beta_interaction)
	
	assign(paste("resultsinterSE_",j,sep=""),se_interaction)

}





# Vectors of Results 	
coef <- 0
coef_interaction<-0
se <- 0
se2 <-0
se_interaction<-0
se_interaction_2<-0



for (j in 1:m){
	coef <- coef + get(paste("resultsc_",j,sep=""))
	coef_interaction<- coef_interaction + get(paste("resultsinteraction_",j,sep=""))
	
	se <- se + get(paste("resultss_",j,sep="")) 
	se2 <- se2 + (get(paste("resultss_",j,sep="")))^2
	
	se_interaction<- se_interaction + get(paste("resultsinterSE_",j,sep="")) 
	se_interaction_2<- se_interaction_2 + (get(paste("resultsinterSE_",j,sep="")))^2	
}


coef <- coef/m
coef_interaction<-coef_interaction/m
se <- se/m
se2 <- se2/m
se_interaction_2 <-se_interaction_2/m
se3 <- 0
se4 <- 0




for (j in 1:m){
	se3 <- se3 + (get(paste("resultsc_",j,sep="")) - coef)^2
	se4 <- se4 + (get(paste("resultsinteraction_",j,sep=""))-coef_interaction)^2

}


se3 <- se3/(m-1)
se4 <- se4/(m-1)
final_se <- sqrt(se2 + se3*(1+(1/m)))
t <- coef/final_se

final_se_interaction<-sqrt(se_interaction_2 + se4*(1+(1/m)))
t_interaction<-coef_interaction/final_se_interaction




# Output Results
cbind(coef,final_se,t)

# Marginal Effects 
cbind(coef_interaction,final_se_interaction,t_interaction)




##### PLOT#####

vkrange=c(1,2,3,4,5)


plot(vkrange, coef_interaction,"l",lwd=2,ylim=c(-.35,.35),xlab="Corporatism",ylab="Marg. Effect of Left Gov",main="GovEm")

lines(lowess(vkrange,coef_interaction + 1.96*final_se_interaction),col="red",lwd=1,lty=3)

lines(lowess(vkrange,coef_interaction - 1.96*final_se_interaction),col="red",lwd=1,lty=3)

abline(h=0,col="blue",lwd=.5)




